/*
 * index2.c
 *
 *  Created on: Jan 29, 2014
 *      Author: jiang
 */


/*
 * index.c
 *
 *  Created on: Nov 8, 2013
 *      Author: jiang
 */



#include <time.h>
#include <inttypes.h>
#include "global.h"
#include "seqio.h"
#include "bntseq.h"
#include "lookup.h"
#include "util.h"

//DEBUG ONLY///////////////////////////
//uint32_t gloablCount = 0;
//uint64_t NinReadsCount = 0;
//uint32_t gloabl3KmerCount = 0;
//uint32_t unableMapReadsCount = 0;
//uint64_t debug_totalKmer =0;
//////////////////////////////////////

/**
 * Read two bits from one byte. The left most index is 0.
 * The Possible values of pos are: 0, 2, 4, 6
 */
static uint32_t read2bits(uint8_t byte, uint8_t pos) {
	uint32_t twoBits;
	twoBits = ((byte) >> (6 - pos)) & ((1 << (2)) - 1);

	return twoBits;
}

int dmap_index_ref2(const char*refFile, opt_t *op) {
	kseq_t *seq;
	uint8_t *pac;
	bntseq_t *bns;

	int64_t l_pac;
	char *prefix = refFile;
	int ret, forward_only = 0 ,testIter, num , startPos;
	uint8_t *reverpac;
	int forward = 0, backward =1;
	int arrayTypeTable =1;
	char *tableFn, *indicatorArrayFn;
	clock_t t;

	uint8_t *kmerStrandIndicator; // bit set 1 means that distinct kmer is distinct  on whole + and - strand, but located in - strand
	table_t *kmer_position_table;


    t = clock();
    fprintf(stderr,"[dmap_index] indexing");
	//initialize
	if (initTableArray1(&kmer_position_table,&kmerStrandIndicator) != 1) exit(1);

	//make reference sequence into compact representation, 2 bit per base
	gzFile fp = xzopen(refFile, "r");

	//l_pac, length of sequence. 4 times the compact representation
	l_pac = bns_fasta2bntseq(fp, prefix, forward_only);

	bns = bns_restore(prefix);
	pac = calloc(bns->l_pac / 4 + 1, 1);
	if (pac == 0) {
		fprintf(stderr, "[dmap_index] insufficient memory!\n");
		exit(1);
	}
	//index reference, core
	{
		//Get compact presentation of sequence into memory
		{
			fread(pac, 1, bns->l_pac / 4 + 1, bns->fp_pac);

			/*
			//TEST reverse complement correctness
			fprintf(stderr,"bns->l_pac / 4 + 1: %d \n",bns->l_pac / 4 + 1);
			for ( testIter = 0; testIter < bns->l_pac / 4 + 1; testIter++)
				fprintf(stderr, "TEST  pac[%d]: %x   \n", testIter,pac[testIter]);
			*/
			fprintf(stderr, "[dmap_index] length of pac %llu \n", bns->l_pac);

			//Build for + strand
			startPos = 0;
			buildTable(kmer_position_table,kmerStrandIndicator,pac,(bns->l_pac)/2, startPos, forward);


			//Build for - strand
			num = ((bns->l_pac)/2)%4;
			if (num == 0)//Good
				startPos = 0;
			else
				startPos = 4 -num;

			//fprintf(stderr,"startPos %d \n",startPos);
			reverpac = pac+((bns->l_pac)/2)/4 + 1 ; //That will ommit startPos number of character of the - strand

			//TEST reverse complement correctness
			/*
			fprintf(stderr,"bns->l_pac / 4 + 1: %d \n",bns->l_pac / 4 + 1);
			for ( testIter = 0; testIter < bns->l_pac / 4 + 1; testIter++)
						fprintf(stderr, "TEST  pac[%d]: %x   \n", testIter,pac[testIter]);

			for ( testIter = 0; testIter < (bns->l_pac/4 +1) - ((bns->l_pac)/2)/4 -1 ; testIter++)
				fprintf(stderr, "TEST  reverpac[%d]: %x   \n",testIter ,reverpac[testIter]);
			*/

			buildTable(kmer_position_table,kmerStrandIndicator,reverpac,((bns->l_pac)/2 - startPos) ,startPos,1);

			t = clock();
			tableFn = (char*)calloc(strlen(refFile) + 10, 1);
			indicatorArrayFn = (char*)calloc(strlen(refFile) + 10, 1);
			strcpy(tableFn , refFile);strcat(tableFn,".tb");
			strcpy(indicatorArrayFn , refFile);strcat(indicatorArrayFn,".idc");

			fprintf(stderr, "[dmap index] write index results to file...\n" );
			if (table_dump(tableFn,kmer_position_table, arrayTypeTable) == 0 && indicatorArray_dump(indicatorArrayFn,kmerStrandIndicator) ==0)
			{
				fprintf(stderr, "[dmap index]  write success! %.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
			}
			else
			{
				fprintf(stderr,"[dmap index] write index file failed \n");
				exit(1);
			}
		}//compact
	}
	return 0;
}

void buildTable2(table_t *kmer_position_table,uint8_t *kmerStrandIndicator,uint8_t *pac, int64_t l_pac, int64_t startPos, int foward_backward)
{
	int i, j;
	uint32_t pos;
	uint8_t test = 0x36;
	uint32_t first, oldKmer, nextKmer;
	uint32_t mask = (1 << ((KMER_LEN -1) * 2)) - 1;

	if ((KMER_LEN -1) == 16 ) mask = 0xffffffff;

	//Build kmer-position mapping. kmer to number, then use that number as array index.

	//first kmer
	first = 0;
	for (i = 0; i < ((KMER_LEN -1) * 2) / 8; i++) {
		first = (first << 8);
		first += (uint32_t) pac[i];
	}
	for (j = 0; j < ((KMER_LEN -1) * 2) % 8; j = j + 2) {
		first = first << 2;
		first += read2bits(pac[i], j);
	}
	if (foward_backward == 0 )
	{
		fprintf(stderr,"[index + strand ]");
		putValueArray1(kmer_position_table, first, POSITION_ZERO);
		//fprintf(stderr,"forward  first\n");
	}
	else if (foward_backward == 1)
	{
		fprintf(stderr,"[index - strand ]");
		if (isEmptyArray1(kmer_position_table,first))//State 0
		{
			putValueArray1(kmer_position_table, first,POSITION_ZERO);
			setKmerStrandIndicatorBit(first,kmerStrandIndicator);
		}
		else if (isNotuniqueArray1(kmer_position_table,first))//State NOT_UNIQUE
		{
		} else if (isUniqueArray1(kmer_position_table,first)) //State 1
			putValueArray1(kmer_position_table, first,NOT_UNIQUE);
		else
		{
			fprintf(stderr,"[Index Ref] ERRO! \n");
			exit(1);
		}
	}
	else
	{
		fprintf(stderr,"[Index Ref] ERRO! \n");
		exit(1);
	}



	//now pac[i] hold byte that we are going to operate, j is starting position(0, 2, 4, 6) of next 2 bits
	oldKmer = first;
	//DEBUG ONLY
	fprintf(stderr, "l_pac - (KMER_LEN -1):  %"PRId64"  \n",
			l_pac - (KMER_LEN -1));
	//DEBUG ONLY

	//position is ZERO based
	for (pos = 1; pos < (uint32_t) (l_pac - (KMER_LEN -1)); pos++) //read  bns->l_pac -1 times, pos = 1  means kmer 1, first kmer  0 start with position 0
	{
		if (pos == 1000000)
			fprintf(stderr, "[dmap_index] %llu characters processed \n", pos);
		if (pos % 50000000 == 0)
			fprintf(stderr, "[dmap_index] %llu characters processed \n", pos);
		oldKmer <<= 2;
		nextKmer = (oldKmer + read2bits(pac[i], j));
		nextKmer &= mask;

		j += 2;
		if (j == 8) {
			i = i + 1;
			j = 0;
		}

		//FIXME  Cross + strand and - strand, do not count

		//if((pos > bns->l_pac/2 - (KMER_LEN -1))  &&  (pos < bns->l_pac/2) ) //Cross + strand and - strand, do not count
		//{
		//	oldKmer = nextKmer;
			//fprintf(stderr,"%llu do not count!\n",pos);
		//	continue;

		//}

		/*Do some processing on kmer*/
		if (isEmptyArray1(kmer_position_table,nextKmer))//State 0
		{
			putValueArray1(kmer_position_table, nextKmer,pos+startPos);
			oldKmer = nextKmer;
			if (foward_backward == 1) // In - strand state 0  ---> state 1
			{
				setKmerStrandIndicatorBit(nextKmer,kmerStrandIndicator);
				continue;
			}
			else
				continue;
		} else if (isNotuniqueArray1(kmer_position_table,nextKmer))//State NOT_UNIQUE
		{
			oldKmer = nextKmer;
			continue;
		} else if (isUniqueArray1(kmer_position_table,nextKmer)) //State 1
		{
			putValueArray1(kmer_position_table, nextKmer,NOT_UNIQUE);
			oldKmer = nextKmer;
			continue;
		}
		else
		{
			fprintf(stderr,"[Index Ref] ERRO! \n");
			exit(1);
		}
	}


}

void static newTable(table_t *kmer_position_table, uint8_t *kmerStrandIndicator,uint8_t *pac, int64_t l_pac)
{
	fprintf(stderr, "compute dark region ... \n");
	//for (index = 0; index < KMER_TABLE_SIZE; index++) {

	//}

}



